A one-dimensional dipole lattice model for water in narrow nanopores 
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We present a recently developed one-dimensional dipole lattice model that accurately captures 
the key properties of water in narrow nanopores. For this model, we derive three equivalent repre- 
sentations of the Hamiltonian that together yield a transparent physical picture of the energetics of 
the water chain and permit efhcient computer simulations. In the charge representation, the Hamil- 
tonian consists of nearest-neighbor interactions and Coulomb-like interactions of effective charges 
at the ends of dipole ordered segments. Approximations based on the charge picture shed light on 
the influence of the Coulomb-like interactions on the structure of nanopore water. We use Monte 
Carlo simulations to study the system behavior of the full Hamiltonian and its approximations as a 
function of chemical potential and system size and investigate the bimodal character of the density 
distribution occurring at small system sizes. 
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I. INTRODUCTION 

Confinement of water into molecularly narrow pores 
dramatically influences its structure and dynamics. For 
instance, the water flow through single wall carbon nan- 
otubes is orders of magnitude larger than that expected 
from continuum hydrodynamics with water fluxes 
comparable with those of transmembrane proteins • In 
carbon nanotubcs with sub-nanometer diameters, water 
forms wires that provide ideal routes for proton conduc- 
tion d, H, [1]. Due to these properties, water-filled car- 
bon nanotubes are promising building blocks for future 
high flux membranes [H, 0, UH, 0, desahnation systems 
[IB 113, mi , and nanofluidic devices , as well as mod- 
els for biological water and ion transport [l^, [3 ■ 

Several recent computer simulation studies have fo- 
cused on the structure and energetics as well as the trans- 
port properties of water in the interior of narrow carbon 
nanotubes [1, 0, S [3 • These simulations indicate that in 
a water bath at ambient conditions such nanotube chan- 
nels fill with water despite the hydrophobic nature of 
their walls a prediction which has been confirmed 
by experiment p^ . [iGj . For sufficiently small nanotube 
diameters, the smooth inner tube wall forces the water 
molecules to stay near the tube axis such that they form 
a hydrogen bonded and orientationally ordered single file 
chain. Remarkably, these chains remain ordered to al- 
most macroscopic length scales of > 0.1 mm due to the 
relatively high energy of orientational defects [l3|- For 
larger system sizes these defects destroy the order and 
no true order/disorder phase transitions occurs in the 
thermodynamic limit as required by theory (Tsl . fioj . 

Due to the reduced mobility of the water molecules in 
the one-dimensional confinement, single file water chains 
in narrow pores can be modelled using lattice models 
with discrete degrees of freedom. Such simplified models 
capture the essential physics of diverse phenomena rang- 
ing from tube filling to the protonic conduction and water 
diffusion B HO, [H, HI. Recently, we have introduced 
such a one-dimensional lattice model, in which water 



molecules are represented 0jS 0jS point dipoles oriented ei- 
ther parallel or orthogonal to the tube axis [l,[l3| • In con- 
trast to other lattice models, our dipole model quantita- 
tively reproduces the structure of quasi one-dimensional 
water in the tube interior including the formation of de- 
fects and their interactions. Moreover, simulations of this 
model are computationally inexpensive such that studies 
of large systems are possible that would not be feasible 
otherwise. 

Here, we present a detailed derivation of our lattice 
model and its various different, mathematically equiva- 
lent, representations. In this model, dipoles are arranged 
on a regular Id-lattice and interact via 1/r^ dipole-dipole 
interactions. This dipole picture can be simplified by 
grouping domains with equal orientation into segments. 
In the resulting segment picture, the total energy of the 
system is written as a sum of the internal energies of the 
segments and their interactions, which are of the dipole- 
dipole type. As we shall see, this segment picture is es- 
pecially useful for the formulation of Monte Carlo moves 
that satisfies the configurational constraints dictated by 
the model. Resummation of the internal energy of the 
ordered segments finally leads to the charge picture, in 
which the total energy is expressed as a sum of Coulomb- 
like interactions of effective charges placed at the end- 
points of the ordered segments. In this physically ap- 
pealing picture the Coulomb- like interactions account for 
all effective interactions of defects, chain ends, and pro- 
tons. Because of its reduced computational complexity, 
the charge representation permits simulations of tubes of 
macroscopic length and investigations of the approach of 
the thermodynamic limit. The charge picture lends itself 
to approximations in which the long-range interactions 
of the effective charges are neglected and which allow the 
investigation of the role of the Coulomb-like interactions. 
In these approximations the filling of the tube is repro- 
duced correctly, but the defect number at the filling tran- 
sition is not. For small system sizes the particle number 
distribution function is bimodal with peaks at low and 
high densities and is captured nicely by the approxima- 



2 



tions. However, we find that neglecting Coulomb- like in- 
teractions qualitatively alters the form of the low density 
peak. 

The remainder of this paper is organized as follows. In 
Sec. HI] we develop the dipole model and derive the two 
other equivalent formulations of the Hamiltonian. One of 
these, the charge picture, is then used to obtain approxi- 
mations of the Hamiltonian. In Sec. Illll we introduce the 
simulation methods and in Sec. IIVI we discuss the pa- 
rameterization. The filling/emptying transition and the 
bistability of the particle number distribution are pre- 
sented and discussed in the context of the approxima- 
tions in Sec. [V] The paper concludes with a discussion 
in Sec. El 
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FIG. 1: ID water wires in nanopores. (a) Chain configura- 
tion with a D-defect and an L-defect, and the corresponding 
lattice model in the (b) dipole representation, (c) the segment 
representation, and (d) the efi'ective charge representation. In 
the lattice model, sites are marked by filled circles, dipoles 
are represented by short arrows, segments by long arrows, 
and effective charges by circles with their sign at the center, 
(a) The D-defect molecule accepts two hydrogen bonds from 
the two neighboring water molecules; in contrast, the L-defect 
molecule donates two hydrogen bonds, (b) The next neigh- 
bor dipoles of the defect sites point to the defect site for the 
D-defect and away from it for the L-defect. (c) The configu- 
rational energy of the water wire due to the dipoles is deter- 
mined by the length, orientation, and distance of segments, 
(d) In the effective charge representation, the segments are 
replaced by charges at the their ends with signs according 
to their orientations. As a consequence, defects are formed 
by pairs of equal charges which are positive for the D- and 
negative for the L defect. 

Water in nanopores forms hydrogen-bonded chains of 
water molecules. Such a chain is dipole ordered if all 
molecules accept a hydrogen bond from the neighbor- 
ing molecule on the left and donate a hydrogen bond 
to a molecule on the right (or the other way round). 



This order is destroyed by orientational defects. Note, 
however, that hydrogen-bond defects within a contigu- 
ous chain preserve the total number of hydrogen bonds. 
Figure[lja) shows a chain of water molecules that consists 
of three ordered segments connected by defect molecules. 
The segments consist of two molecules each and are orien- 
tationally ordered. The D-defect connects two segments 
pointing towards each other and an L-defect connects two 
segments pointing away from each other. In contrast to 
molecules within ordered segments which donate a sin- 
gle hydrogen bond and accept a single hydrogen bond, 
the D-defect molecule accepts two hydrogen bonds with- 
out donating any and the L-defect donates two hydro- 
gen bonds without accepting any. This also means, that 
defects cannot be located at chain ends. Another con- 
figurational constraint is that defects can not be located 
next to each other, because the corresponding molecular 
configurations are unstable and lead to immediate recom- 
bination of the defects. Note, that the defect structures 
shown in Fig. [TJa) are the most typical configurations, 
but others can occur as well. 

The free energetics of such a chain of water molecules 
are captured with great accuracy by a one-dimensional 
dipole lattice model. Molecular simulations show that 
due to the hydrogen bonds the water molecules are on 
average located on the sites of a one-dimensional lattice 
and that the long-range interaction of water molecules in 
an ordered chain is given by the dipole-dipole interaction. 
The magnitude of the dipoles is given by the average of 
the component of the dipole moment along the tube axis 
of a water molecule in an ordered chain. As a conse- 
quence, in the dipole model an ordered segment consist 
of equally oriented dipoles parallel to the tube axis, lo- 
cated on the site of a one-dimensional lattice [Fig.[IJb)]. 
The dipole moments of defect molecules are on average 
perpendicular to the tube axis and we only include their 
next-neighbor interactions. 

On this basis, we can formulate an effective Hamil- 
tonian with a reduced number of degrees of freedom. 
This Hamiltonian describes the free energetics of arbi- 
trarily filled tubes that in general consist of ordered or 
disordered hydrogen-bonded chains of water molecules 
(or fragments), with gaps between them. 

Let us assume that our lattice has N sites. The lattice 
spacing is a and the dipole moment p. A dipole located on 
site v has a direction cr^ = 1 if it points "up" , and ct^ = 
— 1 if it points "down" the tube axis. The interaction 
potential of two dipoles located on sites f and ^ separated 



by a distance di/^ 



fj,\a ^ is given by 



(1) 



If a site v carries a defect or if it is empty, we assign 
this site cr^ = 0. In these cases Eq. ((!]) remains valid as 

= 0. Here, we introduced the energy e = j§^^, 
where eg is the dielectric constant, that sets the scale for 
the dipole-dipole interaction. In the following, we use 
reduced units, i.e., e as unit for the energy, a as unit for 
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the length, and p as unit for the dipolc moment. 



within the segment, i.e. , 



li-i li 



A. Dipole picture 

Molecular simulations show that the interaction en- 
ergy of next neighbor molecules, the so-called contact 
energy, is different from the dipole-dipole interaction en- 
ergy of next neighbors, which is given by 4>i,.v+i = ~1 for 
(^vO'w+i = 1. It will be useful to include the next neighbor 
dipole-dipole interaction in the Hamiltonian and correct 
for it to get the right contact energies. Let the configu- 
ration consist of n water molecules (occupied sites), Uc 
hydrogen bonded chains of molecules (fragments), and 
rid defects. We add the contact energy for each of the 
(n — nc) hydrogen bonds and subtract the next neighbor 
dipole interaction for the {n — Uc — 2nd) next neighbor 
pairs of parallel dipoles. This leads to the following ef- 
fective Hamiltonian for water in nanopores 



H 



N-l N 

E E 



(tiuti + {n-nc){l + Ec)-2nd + ncSc , (2) 



where the double sum extends over all pairs of sites. 5c 
is an entropic contribution that accounts for the different 
contributions to the phase space volume of molecules at 
the chain ends and at defect sites compared to molecules 
within the chain. These different contributions are re- 
lated to the number of dangling OH bonds (see Sec. IIV|) . 
As we will see, this Hamiltonian is just one of three equiv- 
alent descriptions of the system. We will refer to this way 
of calculating the Hamiltonian as the dipole picture [see 
Fig. [TJb)], as we sum over all dipole-dipole interactions. 



B. Segment picture 

We can also formulate the Hamiltonian of Eq. ([2|) in 
terms of the Ug = Uc + rid segments, numbered from left 
to right, leading to the so-called segment picture [see Fig. 
He)]. 

For this purpose, we need the dipolar internal energy 
of a segment, which stems from the interaction of the 
dipoles within a segment, and the dipolar interaction en- 
ergy of two segments, which stems from the interaction of 
dipoles belonging to two different segments. For brevity 
we will drop the term "dipolar" and only speak of inter- 
nal and interaction energies of segments in the following. 

The beginning of the segment with index i is given by 
its coordinate Xi and its length is denoted by li. The 
coordinates of the dipoles of segment i are given by Xi + 
n — 1/2 with 1 < n < li. As all dipoles within a segment 
i have the same direction, i.e., ai, = Si for all dipoles 
within the segment i, we assign each segment a direction 
Si ~ ±1. The internal energy E{li) of a segment i is 
given by the sum over all pair interactions of dipoles v 



m)--T. E 



(3) 



The interaction energy I{li,lj, SiSj, Alij) = lij of seg- 
ments i and j, with i < j, is defined as 
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(4) 



with the size of the gap between the two segments given 
by Alij = Xj - {xi + k). 

In the segment picture we obtain for the Hamiltonian 
of a single chain 



H = Ei?(;.)+ E E ^^^ + 

+ {n - nc){l + Ec) - 2nd + UcSc 



(5) 



The calculation of the energy according to Eq. ([5]) in 
the segment picture can be simplified considerably by de- 
riving an explicit functional form for the internal energy 
and expressing the interaction energies in terms of the 
internal energy. Firstly, we can avoid the double sum in 
the calculation of the internal energy of Eq. ^ by count- 
ing all pairs of dipoles separated by a certain distance. In 
a segment of length I, the interaction potential of dipoles 
separated by a distance j with 1 < j < I appears {I — j) 
times. This leads to 
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i~i I- 

Er^-^E^- 
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(6) 
(7) 



for Eq. ©. 

For Z — > cx), the two sums in Eq. ([7]) the two sums can 
be expressed in terms of Riemann's zeta function [2^ 



dm) = E-^' 



(8) 



Accordingly, for long segments we can approximate the 
internal energy E{1) as 



i?z»i(0 = c(2)-;c(3) 



(9) 



The difference, $(/) = Ei^i{l)—E{l) , between the above 
approximation and the exact internal energy is given by 
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m = Ei-^Ei 



and can be rewritten as 



m = *'(o + 2*"(o 



(10) 
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FIG. 2: Calculation of the interaction energy. The interac- 
tion energy / of configuration (e), consisting of two chains of 
length h and Ij, separated by a distance, Ahj, of four sites, 
can be calculated by subtracting the internal energies E of 
configurations (b) and (c) from the internal energy of (a), 
and adding the internal energy of the configuration (d) . 



using the polygamma function [22 



(I+J) 



7n-\-l 



(12) 



Hence, the internal energy of a segment of length I can 
be written as the sum of a linear part, Eq. and a 
non-linear part, Eq. pT|) . leading to the exact expression 



Eil) = c{2) - Kis) - m 



(13) 



The next step towards a simpler energy calculation in 
the segment picture is to express the interaction energy 
of two segments i and j in terms of internal energies (see 
Fig. [2) . As the directions of the segments only determine 
the sign of their interaction energy, we assume, for the 
moment, that the two segments of length li and Ij have 
the same direction, Fig. [21(e). 

In the following we use the simple fact that the internal 
energy E{1) of a chain of length / can be calculated from 
its parts of length /' and I" , with I = I' + I", as 



E{l) = E{l') + Eil") + 1(1', I", 1,0) 



(14) 



For brevity we drop Sij and Alij as arguments of 
I{li,lj,Sij,Aij) = here. Using the equation 

above the internal energy of a segment with length I = 
h + Ij + Alij [Fig. [H (a)] can be written as 



E{1) = Eih) + E{lj) + E{Ak,) + 

+ I{kJ,) + Iik,Ak,) + I{Akj,lj) 



(15) 



To get the desired interaction energy, Ij) on the right 
hand side of Eq. ([TSj) , we subtract the internal energies of 
a segment of length li + Alij and of a segment of length 
Alij + Ij [Fig. [11(b) and (c)] . These two energies, given 
by 



Eik + 
E(Akj 



Al 



h) 



E{k) + E{Akj)+I{k,Akj) (16) 
E{Akj) + E{lj) + I{Ak,,l,) (17) 



include the internal energies of the two segments i and j, 
and the interaction energies of the segment i and j with 
the segment of length Alij separating them. They also 
include twice the internal energy E{Alij) of the segment 
in the middle. To correct for this double subtraction, we 
have to add this energy once. Fig. [2l (d). As the linear 
terms of Eq. (fTS]) cancel we obtain for the interaction 
energy 

I,j = -S^S■j [$(?, + Ij + AUj) + $(A/y) + 

-•^{U + Akj)-^{lj+AUj)] . (18) 

Note that this expression for the interaction energy is 
also valid for continuous gap distances Alij > 0. 

This result can also be obtained by simply inserting 
the expressions for the interaction energies I{li, Alij) and 
I {Alij, I -i) following from Eqs. HH) and ([T7|). 

I{k,Akj) = E{k + Ak,)-E{li)-E{Akj) (19) 
I{Akj,lj) = E{Akj +lj)- E{Akj) - E{lj) (20) 

to Eq. P^. 

Using Eqs. ([T51) and ([T51) we can write the total energy 
in the segment picture, Eq. ([5]), as 

Jls — 1 "s 

i=l j=i+l 

-$(?. + Ah,) - $(/, + Akj)] - -^ih) + 

i=l 

+ndCd + ricCc + nc . (21) 

This Hamiltonian is linear in the number of defects rid, 
the number of fragments ric, and the number of occu- 
pied sites n. The corresponding coefficients Cd, ct, and c, 
which depend on the contact energy E^ and the entropic 
contribution Sc, are given by 



Cd = C(2) + C(3)-2 , 
Cc = Ci2)-l-E, + Sc 
c = 1 + E,- C(3) . 

C. Charge picture 



(22) 
(23) 
(24) 



We can use this expression for the Hamiltonian in the 
segment picture to derive another equivalent description 
of the system, i.e. the so-called charge picture [see Fig. 
[IJd)] . The first few terms of the series expansion of the 
non-linear part of the internal energy. 
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O 



, (25) 



show that for large chain lengths / the non-linear part is 
proportional to l/l, i.e., Coulombic. Thus, we replace all 
segments by charges at their beginnings and their ends 
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according to their orientations. As each segment carries 
two charges, the indices of the charges are given by 2z — 1 
for the charge at the beginning and 2i for the charge at 
the ends of segment i. The charge at the beginning of 
segment i is (721-1 = and the charge at the end is 
q2i ~ Si. The coordinates of these charges are Z2i-i = Xi 
and Z2i = Xi + li. The interaction potential of two charges 
Qm and q„ separated by a distance Zmn = \zn — Zm\ is 
given by ^dzmn) = qraQn^izmn)- This interaction is 
Coulombic in character for large distances and we can 
rewrite this interaction as 

e^cizmn) « e- = . (26) 

with the magnitude of the charges then given by \Qm\ = 
p/a. Figure [3] shows a comparison of this Coulomb-like 
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FIG. 3: Comparison of the Coulomb interaction with the non- 
linear part of the internal energy. While the Coulomb inter- 
action diverges for r ^ 0, "l'(O) — C(2). 

interaction, $(r). with the Coulomb interaction l/r. The 
main difference is that for a distance r — > the Coulomb- 
like interaction converges to a finite value, ^'(0) = 
whereas the Coulomb-interaction diverges to — c». Note 
that for all possible distances, Zmn > 1, the difference is 
small. 

Rearranging the Hamiltonian in the segment picture 
given by Eq. (|2ip we obtain the Hamiltonian in the charge 
picture as 

2na-l 2na 

-ff = ^ ^ grnqn^iZmn) + 
m— 1 n—m+l 

+ndCd + ricCc + no . (27) 

The sum in Eq. ([27|) includes the interaction of the 
charges belonging to the same defect, given by $(1) = 
C(2) — C(3), which occurs times. Introducing the de- 
fect excitation energy, i.e., the energy (the free energy in 
the molecular model) needed to introduce a single defect 
in an infinitely long chain, 

E-D=Cd- 0(1) = 2C(2) - 2 (28) 



we can write the Hamiltonian in the charge picture as 

2na-l 2ris ' 

= ^ ^ qmqn^iZrnn) + 
rn—1 n— m+1 

+n^E'D + jicCc + nc . (29) 

where the prime indicates that the sum does not include 
the interaction of two charges belonging to the same de- 
fect with each other. The coefficient Cc corresponds to 
the energy needed to break an infinitely long chain, and 
move the resulting fragments infinitely far apart from 
each other. The coefficient c is the energy required to add 
a single dipole to an infinitely long chain. The Hamilto- 
nian in the charge picture depends on the charge posi- 
tions, Zj, the number of particles n, the number of chains 
ric, and the number of defects 77.d. We note that these 
quantities are not sufficient to specify a configuration of 
the dipole model unambiguously. 

The Hamiltonian in the charge picture [Eq. ((29|) ] high- 
tlights the Coulomb-like effective interactions of defects 
and chain ends. L- and D- defects attract each other 
Coulombically whereas defects of the same kind repel 
each other. Defects next to a chain end are always at- 
tracted by the charge at this end. 

D. Approximate representations 

If the distances between charges are large we should be 
able to neglect their Coulomb-like interactions that decay 
as l/r. Then, only the linear terms of Eq. (|29p remain. 
The Hamiltonian in this simplest approximation, which 
we refer to as the no-charge- approximation (NCA), is 
given by 

Tio = nc-\- ricCc + ri^Ex) ■ (30) 

For small distances of effective charges the Coulomb- 
interaction is strong and cannot be neglected. This is the 
case for short chains and segments which carry charges at 
their ends. In particular, chains and segments of length 
one are represented by two charges of opposite sign which 
are separated by one lattice constant only. At the next 
higher level of approximation, we add the interaction en- 
ergy — <&(!) = — C(2) + C(3) between the charge pairs 
associated with each of the ni segments of length one. 
The resulting Hamiltonian is 

Hi =nc + ncCe + ndSi, + ni[-C(2)-f C(3)] . (31) 

We will refer to this Hamiltonian as the singlet-charge- 
approximation (SCA) since we include the interaction of 
charge pairs associated with single dipoles. 

One main feature of these approximations compared 
to the full Hamiltonian is that the approximations can 
be written as spin models (with three spin states for the 
NCA and four spin states for the SCA) with nearest- 
neighbor interactions only. 
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The NCA and SCA do not depend on the length, posi- 
tions, or orientations of the segments. States specified by 
the particle number n and no = {7ic,7id} for the NCA, 
and by n and ni = {7ic,nd,ni} for the SCA are thus 
degenerate. Calculating the degeneracy of these states 
requires to count the number of states for these approx- 
imations as a function of these variables. For the NCA 
the number of states, To{N,n,nQ), is given by 

with the number of segments, Ug ~ nc + n^, being a 
function of the defect number and the chain number. For 
the SCA the number of states, ri(A^, n,ni), is given by 



Fi = 2" 



n - 2ns + He - l\ / ns\ - I 
ris - rii - 1 / \nij \nc - 1 
N + l-n 
rir 



(33) 



for n > 0. Uris ^ nc = n then Fq = Fi = 2"(^+|""). 
For a derivation of these equations and see App. |^ 



the last at 72 = ji + ^ ^ 1- Using the polygamma function 
we obtain 



(35) 



00 



00 



- 'T.^'"-^ E r'- (36) 

= s[^-^'{ji) + ^'{:n + i)] ■ (37) 

for the interaction energy of the proton with the segment. 
Rewriting this interaction energy as a function of the 
distance x = | ji — i | -f 1/2 of the beginning of the segment 
from the proton we get 

W{x,l) ^ s[^^'{x + 1/2) + ^'{x + 1 + 1/2)] = 

= s[^{x) + ^{x + l)] (38) 

with a Coulomb-like interaction given by 
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E. Proton defects 

Up to now we only considered chains of intact water 
molecules. If an excess proton is introduced into the sys- 
tem, it forms a hydronium ion consisting of one oxygen 
atom and three hydrogen atoms. In a single-file chain, 
this hydronium ion donates two hydrogen bonds without 
accepting any. Thus, structurally the hydronium ion cor- 
responds to an L-defect with an additional proton Q . In 
the charge picture, the effective interaction of this proto- 
nated L-defect is given by the sum of the effective inter- 
action of an L-defect and the effective interaction of the 
proton with defects and chain ends. 

In Ref. 3 an approximate expression for this effective 
interaction was derived. For completeness we show next 
how the effective interaction of an additional proton lo- 
cated on an L-defect is included in the charge picture. In 
the dipole model the interaction of a proton located on 
an L-defect at site i and a dipole with direction aj — ±1 
at site j is given by 



^ 1 ajpejj - i) 



e'^^ (34) 



where e is the elementary charge and s' = eae/{2p) the 
unit of the energy for the rest of this subsection. 

Using the ideas of the derivation of the segment and 
the charge picture, we calculate the interaction energy 
of a proton and an ordered segment with / dipoles of 
relative orientation s = (Tj{j — i)/\j — i\ with respect 
to the position of the proton, i.e., s = 1 (s = —1) for 
dipoles pointing away (towards) the proton. Let us say 
the proton is located in the origin at the site with index 
zero, the first dipole of the segment is located at ji, and 



Thus, the interaction energy of a proton with the dipoles 
of a segment is equal to a Coulomb-like interaction of 
the proton with charges at the beginning and the end of 
the segment. For a protonated L-defect in an infinitely 
long chain we obtain a proton energy of Ep = —2^(2), 
stemming from the two charges next to the defect. 

For large distances, the interaction energy of a hydro- 
nium ion with a defect can be written as 



«'pl(2) 



1 tQ 



(e-Q) 



(40) 



where Q = 2p/a is the magnitude of the total effective 
charge of a defect. The plus sign is valid for the inter- 
action with a D-defect and the minus sign for that with 
an L-defect. Since the proton charge is larger than Q, 
a hydronium ion is repelled by D-defects and attracted 
to L-defects, in agreement with Ref. Q. Also, the hy- 
dronium interacts repulsively with the endpoints of an 
otherwise ordered water chain. As a consequence, the 
preferred position of an excess proton in an isolated wa- 
ter wire is at the chain center. 



III. MONTE CARLO SIMULATION 

Next, we examine the thermodynamic behavior of wa- 
ter in nanopores in contact with a heat bath and a par- 
ticle reservoir. Thus, we perform canonical and grand- 
canonical Monte Carlo simulations of the dipole model. 

For the Monte Carlo simulations of the dipole model we 
describe configurations in the segment picture in which 
gaps of unoccupied sites are represented as segments with 
s = 0. Thus, a configuration consists of occupied sections 
(segments) and empty sections (gaps). (In the following 
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we use the term "segments" only for segments of chains 
and the term "section" for segments of chains and empty 
gaps.) Between segments one finds either a defect or an 
empty section. The boundary conditions are given by two 
empty gaps at the beginning and the end of the lattice 
(see App. A configuration is unambiguously defined 
by the beginning x^, the length li, and the value of Si for 
each section i. Defects arc located between next neighbor 
sections with s ^ 0. Including the empty end sections, 
there are hq = + 1 sections with s = corresponding 
to gaps. 

With a configuration given by Ci = {{xj,lj, sj} : 1 < 
J < Hg + uq} the canonical partition function can be 
written as 

Z^(/3,n) = ^e-^™. (41) 

{Ci} 

The grand canonical partition function is given by 

oo 

EN{P,z) = l + J2ZN{f3,n)z" . (42) 

n=l 

The fugacity is defined as z = e^^, where /i is the chem- 
ical potential and (3 ~ \/{k-QT) the reciprocal tempera- 
ture with fee being Boltzmann's constant. 

To enhance the sampling, we use non-local Monte 
Carlo moves. Monte Carlo simulations with only local 
moves, in which individual dipoles are flipped, are ineffi- 
cient for large systems. Here, we apply efficient non-local 
trial moves with asymmetric generation probabilities for 
which we have to correct in the acceptance probability. 
These trial moves change the lengths of sections and their 
orientations for the generation and recombination of de- 
fects, for the displacement of defects and chains, and for 
the insertion and deletion of particles. For details see 
App.E 

We also perform Monte Carlo simulations for the NCA 
and SCA for which the canonical partition functions are 
given by 

Z^*)(n) = ^ (43) 

with 1 = for the NCA and i = 1 for the SCA, and 
implicitly depending on the configurations Cj . Using the 
degeneracies given by Eqs. ([5^ and ([55]) . Eq. (|^ can 
be rewritten as 

= E nOe-''^'^"'"') , (44) 

which permits us to formulate new effective "Hamiltoni- 
ans" 

n',^n^-T\nT,{N,n,n,) (45) 
with canonical partition functions given by 

Z«(n)=Ee-^«^("^"-) . (46) 



We can calculate these partition functions either numeri- 
cally by direct summation or perform Monte Carlo simu- 
lations in the space of ng — {ric, rid} for the NCA and of 
ni = {nc,n.d,n.i} for the SCA. The applied trial moves 
simply increase or decrease the values of the variables 
of the Hamiltonian within the limits of the phase space 
given in App. 

To study the system behavior over a broad range of the 
chemical potential we use the Wang-Landau algorithm 
with the particle number as order parameter to find 
a bias function w{n) corresponding to the negative free 
energy as a function of the particle number. We use this 
function for a biased simulation [25| at the fugacity z, 
resulting in a flat histogram of the particle number, with 
the Hamiltonian in the biased system given by 

H' = H -Tw{n)\nz (47) 

The output are samples of the total energy, chain num- 
ber, defect number, number of particles, and the total 
dipole moment, {E^^\ nc \ nl^\ n^^\ D^''^ . By unfolding 
the bias function w{ri) and reweighting [25j we obtain 
estimates for observables that are functions of the above 
quantities (see App. [Cjl . 



IV. PARAMETERIZATION 

The dipole model was parameterized with and vali- 
dated against detailed molecular simulations [l3|. The 
molecular model consists of a (6,6)-type carbon nan- 
otube, filled with up to 100 water molecules, following 
Ref. [2il. We used the TIP3P potential [13 for the water- 
water interactions and a cylindro-symmet ric p otential for 
the tube-water interaction as in Refs. [3, H^. Boundary 
effects are minimized by using tubes that are longer than 
the volume accessible to the water molecules. 

The lattice spacing a, the dipole moment p, and the 
contact energy Ec can be determined in a canonical 
Monte Carlo simulation of a single hydrogen bonded 
chain of water molecules in a nanopore. 

The lattice spacing is given by the average distance of 
neighboring molecules within a chain for which we ob- 
tain a = 2.65 A. The dipole moment is given by the 
average dipole moment of a water molecule in an or- 
dered chain along the tube axis, p = 1.9975 D. This 
determines the energy scale for the dipole-dipole interac- 
tion as £ = 2p^/(47reoa^) = 25.8236 kJ/mol, such that 
(3e « 10.42 at T = 298 K. The contact energy Ec = -20.8 
kJ/mol is determined as the average interaction energy 
of two neighboring water molecules within a chain. 

To perform grand canonical simulations of the dipole 
model we have to determine the effective tube-water in- 
teraction. Since it couples to the particle number we 
can absorb it in the chemical potential. We also need 
the entropic contribution which couples to the number 
of chains. We determine both quantities by comparing 
the transfer free energies of the molecular model and 
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FIG. 4: The transfer free energy as a function of the density 
for different system sizes for the molecular model (solid lines) 
and the dipole model (dashed lines). The excellent agreement 
renders the curves for the different models nearly indistin- 
guishable. Results are for fugacity zq- 



the dipole model and tuning the fugacity and the en- 
tropic contribution to obtain agreement between these 
two models. The transfer free energy, given by jiAin) = 
— ln[P(ri)/P(0)], where P(n) is the particle-occupancy 
distribution function, is the free energy needed to intro- 
duce n water molecules into a tube of length L at the 
reciprocal temperature /?. We obtain for a tube of length 
h = 30a in contact with a heat bath and a particle reser- 
voir at ambient conditions, i.e., at room temperature and 
atmospheric pressure, a fugacity zq = 0.000327 and an 
entropic contribution /35c = —3.96. 

Figure 3] shows a comparison of the transfer free en- 
ergies of the molecular and the dipole model for system 
sizes N = 5, 10, 15, and A'' = 30 as a function of the 
density p = n/N. Since the dipole model is a coarse- 
grained description of the molecular model, all states in 
the molecular model with an occupation number equal 
to or larger than the number of sites contribute to the 
completely filled state, i.e., to P{N). The agreement is 
excellent as the curves for the molecular and the dipole 
model lie practically on top of each other. We observe 
that for small sizes the system shows two minima, one for 
the empty and one for the filled state. Thus, the parti- 
cle number distribution function is bimodal, i.e., has two 
peaks (see Sec. IV Bl) . For growing system size the min- 
imum corresponding to the filled state gets deeper and 
the minimum corresponding to the empty state vanishes. 

By reweighting the particle number distribution func- 
tions corresponding to the transfer free energies shown 
in Fig. [4] we obtain the adsorption isotherms and the rel- 
ative particle number fluctuations as a function of the 
relative fugacity (Figs. [5] and El respectively). The rela- 
tive fugacity is the actual fugacity divided by the fugacity 
of a system in contact with a heat and particle reservoir 
at room temperature and atmospheric pressure. For low 
fugacities the relative fugacity is equal to the relative 
humidity. For larger system sizes the density shown in 



FIG. 5: The particle density as a function of the relative fu- 
gacity. The solid lines are results from molecular simulations 
and the dashed lines are results for the dipole model. 
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FIG. 6: Particle number fluctuations as a function of the 
relative fugacity for different system sizes. The solid lines are 
results from molecular simulations and the dashed lines are 
results for the dipole model. 

Fig. [5] gets steeper at the filling transition, which moves 
to smaller fugacities. 

The relative variance, i.e., the variance of the parti- 
cle number divided by the average particle number, is a 
measure for the fiuctuation of the particle number. For 
macroscopic volumes, the relative variance is related to 
the isothermal compressibility kt via ((n^) — (n)^)/(n) = 
pk^TK.^, where angled brackets denote ensemble aver- 
ages. The relative variance (Fig. [5]) has its maximum 
at the filling transition. The fiuctuations become larger 
with increasing system size. The properties of the fill- 
ing transition for increasing system size are discussed in 
detail in the Sec. |Vl 

The excellent agreement of the dipole model and the 
molecular model for different system sizes supports the 
validity of the entropic contribution. It accounts for the 
different contributions to the phase space volume of wa- 
ter molecules according to the number of their dangling 
OH bonds. In an ordered chain all molecules donate a 
single hydrogen bond except one molecule at one of the 
chain ends. In contrast to all other molecules of the chain 
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PAS 



10° 20° 30° 40° 



A. Drying/Filling Transition 



3.50 2.70 2.15 



TABLE I: Estimates for the entropic contribution. 



which have a single dangling OH bond, this molecule has 
two. Thus, it has more freedom to move and a higher 
contribution to the entropy of the chain than the other 
molecules. If we generate an L-dcfect, which donates two 
hydrogen bonds without accepting any, both molecules 
at the chain ends have two dangling OH, thus conserving 
the number of dangling OH bonds. The situation for the 
D-defect is similar and thus the number of dangling OH 
bonds is conserved for any number of L- and D-defects 
in the chain. Only if a hydrogen bond is broken/formed, 
the number of dangling OH bonds is increased/decreased 
by one. 

Since the entropic contribution accounts for the differ- 
ence in the phase space volume contributions of a dan- 
gling OH bond and one that donates a hydrogen bond, 
we can obtain an estimate from simple geometric con- 
siderations. We assume that an OH bond of a water 
molecule within a segment that donates a hydrogen bond 
is restricted to a spherical cap with an opening angle a. 
For a hydrogen bond of a water molecule at the chain 
end that accepts a single bond, this cap is about half a 
sphere. The ratio of these areas gives us an estimate for 
the difference of the contributions to the entropy of these 
two water molecules. The surface area of a spherical cap 
with opening angle a of a sphere with radius r is given 
hy A = 27r?'^(l — cos a). The difference in the entropies 
between a dangling OH bond and a hydrogen bonded OH 
bond is given by 

The entropic contribution is an energy in our effective 
Hamiltonian and thus equals the negative entropy differ- 
ence, i.e., Sc = —AS*. Table U shows results for different 
opening angles a which are indeed of the order of the 
entropic contribution, pSc = —3.96. 

V. RESULTS 

The charge picture of the dipole lattice model gives 
direct physical insight into the microscopic properties of 
water in nanopores. The Coulomb- like interactions lead 
to an attraction between L- and D-defects and also to 
an attraction between a chain end and the defect next 
to it. So, the Coulomb interaction has clearly an impor- 
tant influence on these microscopic properties of water 
in nanopores. The question arises, what is the influence 
of the Coulomb-like interactions on the overall phase be- 
havior? Or, in other words, which aspects of the system 
behavior are captured by the approximations that neglect 
Coulomb-like interactions (NCA and SCA)? 




FIG. 7: Particle density of the full Hamiltonian (black) and 
the SCA (red). The orange line shows results of NCA for 
iV = 100. 

To clarify the role of the Coulomb-like interactions, we 
compare results derived in the SCA with results of the 
full Hamiltonian for the filling transition. We character- 
ize the system by the particle density, p = n/N, (Fig. [7]), 
the particle fluctuations (relative variance. Fig. [8]), the 
chain density pc = Uc/N, i.e., the number of chains per 
site, (Fig. [5]), and the defect density pd = Ud/N, i.e., the 
number of defects per site (Fig. (TUl), as a function of the 
relative fugacity. The results for system sizes A^=10^, 
10^ 10^, and 10^ sites are obtained by reweighting of 
biased sampling simulation data. Additionally, we show 
results for the SCA for N = 10^° from Monte Carlo sim- 
ulations at the corresponding fugacity values using the 
effective Hamiltonian of Eq. (|45p . 

The adsorption isotherms (i.e., the average particle 
density) in Fig. [7] and the relative fluctuations of the par- 
ticle number in Fig. [5] show that the results for the SCA 
are in excellent agreement with the results for the full 
Hamiltonian for the system sizes studied here. We find 
that both the adsorption isotherms and the relative vari- 
ance are nearly converged to their thermodynamic limits 
for a system size of 10* sites which is supported by the 
results for SCA for N = 10^". 

The adsorption isotherm become steeper with increas- 
ing system size but the slope remains finite even in the 
thermodynamic limit. The relative variance is peaked at 
the filling transition as shown in Fig.[8l Even though the 
peak height initially grows with increasing system size, 
it eventually converges to a finite value. This is in agree- 
ment with the impossibility of a true first-order phase 
transition in one dimension for 1 /t^ interactions [18| . 

Figures [7] and [5] also show results for the NCA for a sys- 
tem size = 100. The adsorption isotherms are in good 
agreement, but the relative variance decays too slowly 
for fugacities below the filling transition compared to the 
results of the full Hamiltonian. The SCA reproduces the 
results of the full Hamiltonian well, as the system at low 
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fugacities consists mainly of chains of length one which 
are correctly described in the SCA. 

Since the adsorption isotherms for the NCA and the 
full Hamiltonian are in good agreement we use this ap- 
proximation to obtain an estimate for the chemical po- 
tential //i/2, where the system is half full. Assuming 
that no defects exists, and that the number of particles 
is much larger than the number of chains we obtain a 
Hamiltonian that only depends on the particle number, 



Ti! = nc — , 



(49) 



where we treat the chemical potential like a magnetic 
field in the Ising model and put it in the Hamiltonian. 
The n particles do not interact with each other but couple 
to the field c — /i. The canonical partition function of 
this ideal lattice gas in an external field is given hy Z — 



(1- 



-/3(c-M)^w 



) . The density is obtained by p - 



-l/ZdZ/d{Pjj,) which gives p = e 



-/3(c- 



V(l^ 



(n) /N : 

-/3(c-At) 



If we demand that the system is half full, i.e. , the density 
is p = 1/2, then the energy of putting a particle into the 
system is equaled by the chemical potential, p = c, which 
leads to 



Pi/2 + C(3) 



(50) 



The corresponding relative fugacity ^1/2/^0 ~ 0.0841 is 
in good agreement with the simulation result ~ 
0.084, corresponding to 8.4% relative humidity. 

This result can also be derived for a lattice gas 
with dipole-dipole interactions and exploiting the isomor- 
phism to the Ising model (sec Add. [P]) . 
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FIG. 8: Relative variance of the particle number of the full 
Hamiltonian (black) and the SCA (red). The orange line 
shows resuhs of NCA for N = 100. 

Observables depending on the particle number are well 
reproduced in the SCA. We obtain insight into the struc- 
ture of the system by looking at the chain and defect 
density with changing chemical potential. We find that 
at the filling transition both the chain density in Fig. [2] 
and the defect density in Fig. [10] are peaked. The peak 
in the chain density is also reproduced by the SCA al- 
though the curves are below the results of the full Hamil- 
tonian. Thus, the peak in the fragment density reflects 




FIG. 9; Chain density of the full Hamiltonian (black) and the 
SCA (red). The blue line shows the density of non-interacting 
particles. 
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FIG. 10: Defect density of the full Hamiltonian (black) and 
the SCA (red). 



the entropic gain through fragmentation for low particle 
densities. 

Figures O and [TO] also show that for the full Hamil- 
tonian the defect density roughly mirrors the fragment 
density. The reason is, that on the one hand a chain with 
a defect can lower its free energy by splitting the chain 
into two at the defect site. In contrast, the generation 
of defects costs less energy at the ends of water chains 
and the number of ends is proportional to the number of 
fragments. The latter is a consequence of the Coulomb- 
like attraction between the charge at the chain end and 
charges forming a defect next to this end. Thus, the re- 
sults of the SCA do not show a peak in the defect density 
and do not reproduce the results of the full Hamiltonian 
at the filling transition even qualitatively. 

For fugacities z below the filling transition, the aver- 
age fragment number decays approximately linearly with 
z towards zero. At low fugacity, the system behaves ide- 
ally. The grand-canonical partition function of ideal par- 
ticles is given hy Z = J2n Cn)i'^^^~^^°)" which leads for 
small fugacities to a density p « 22; exp(— /3S'c), where Sc 
is the entropic correction. The fragment number is then 
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approximately equal to the number of particles, explain- 
ing the observed linear decay for small z. Moreover, with 
fragments too short to carry a defect, the average defect 
number decays even faster. 

The increasing fragmentation while approaching the 
filling transition is also the reason for the faster decay of 
the orientational order for larger system sizes. As mea- 
sure for the order we use the average of the total dipole 
moment squared divided by the system size, (£)^)/A^^ as 
shown in Fig. [TT] for different system sizes as a function 
of the relative fugacity. This order parameter is close to 
unity if the system is orientationally and translationally 
ordered and approaches zero for increasing disorder. Al- 
though the filling transition is steepest for = 10^ the 
order parameter decreases the fastest. This is direct con- 
sequence of the peak in the chain density and the larger 
gaps for this system size which leads to a weaker coupling 
of the chains. 

This conclusion is supported by results for the SCA 
(also shown in Fig. fTTj) which agree nicely with results for 
the full Hamiltonian for small system sizes of Af = 100 
and A^ — 1000 where fragmentation plays a minor role. 
For larger system sizes we observe that the square of the 
total dipole moment is lower for the SCA then for the 
full Hamiltonian. This indicates, that for the full Hamil- 
tonian chains are coupled to their next neighbors via 
Coulomb-likc interaction. Since the SCA is lacking these 
interactions, chains are uncorrclated leading to a lower 
expectation value for the total dipole moment squared. 
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FIG. 12: Bistability map. Density plot of the bistability mea- 
sure M of Eq. (|62p . which is close to one for a bistable system, 
for all investigated system sizes (top) and an enlarged view 
for large systems (bottom). The black solid lines are lines of 
constant density p = 0.01, 0.5 and 0.99 from left to right. 



FIG. 11: The average value of the total dipole moment 
squared divided by the system size squared as a function of the 
relative fugacity for different system sizes for the full Hamil- 
tonian (solid lines) and for the SCA (dashed lines). 



B. Bistability 

It was already observed for small systems at ambient 
conditions that the particle number distribution function 
is bimodal [ij, US ■, showing one peak for the empty and 
one for the full tube. In the following we investigate how 
this bistable behavior changes with system size and chem- 



ical potential and discuss the influence of the Coulomb- 
like interactions. 

Similar to Maibaum and Chandler in Ref. [2^ we start 
the discussion with the Hamiltonian of a one dimensional 
lattice gas in an external field /i, 

7V-1 N 

H = s^s^+i ~ h'^s^ (51) 

i=l i=l 

where the occupation number of site i is given by Si = 0, 1 
and the coupling constant of occupied sites by J > 0. 
This Hamiltonian can be written as 

H = -{J + h)n + Jn^ (52) 
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where n is the total occupation number and the num- 
ber of domains (or chains) of particles. The above Hamil- 
tonian is of the form 

H = ]Cn + Trie (53) 

where /C is the coupling constant and X the interface 
energy. The latter is the energy needed to break a chain 
into two and form two new interfaces between occupied 
and empty sites. In case of the NCA of the dipole model, 
the coupling constant is given by /C = /3c — In z and the 
interface energy by T = (3cc — In 2. 

The particle number distribution of an ideal lattice gas 
(i.e. X = 0) is given by the binomial distribution and 
therefore shows only a single peak. Only with a posi- 
tive interface energy bistability of a partly or completely 
filled system and the empty system is observed. Then 
the particle number distribution function is given by 

P{n) oc e^^'^" J2 "c)e"'^"'^^ (54) 

Tic 

where r(7T,, ric) is the number of states depending on the 
particle number and the number chains, and a functional 
form that depends on the boundary conditions. 

For periodic boundary conditions the number of states 
is given by 

^^^'"'"^ (ric - l) ( ric )^( ric-l ) (ric) ' 

For a constant number of chains, this function is uni- 
modal and symmetric with respect to the location of its 
maximum at n = N/2. Thus the particle number dis- 
tribution function is unimodal if one excludes the empty 
state. Free rather than periodic boundary conditions in- 
troduce a small asymmetry in the number of states, i.e. 

which is not sufficient to change this behavior. 

For a positive interface energy all non-empty states 
are energetically penalized according to their number of 
chains compared to the ideal lattice gas. Thus, the weight 
of the empty state in the partition function increases rela- 
tive to the non-empty states and a second peak for 7i = 
appears (see Fig. [T^ . 

In contrast to this, water in nanopores shows a low 
density peak at densities larger than zero. The reason 
is that the Coulomb-like interaction of charges of oppo- 
site sign at the ends of short ordered chains lowers the 
internal energy of such chains compared to the energy 
they would have with only nearest-neighbor interactions. 
This is already the case for the SCA, where chains of 
length one, i.e., single dipoles for defect free systems, 
are treated separately. In the following, we quantify this 
bistable behavior for the SCA as a function of the system 
size and the fugacity. We calculated the particle number 




n 



FIG. 13: Particle number distributions for different values of 
the interface energy for a system with free boundary condi- 
tions of size A'^ = 100 close to the filling transition {JC = 0). 
For an interface energy 7 > 6 we see a second peak for the 
empty system (n = 0). 

distribution function for system sizes up to iV = 1000 for 
a certain value of the fugacity z. The particle number 
distribution fmiction for the SCA is given by 

P(n) oc z"Z^^'(n) (57) 

with the canonical partition function Z^j^\n) given by 
Eq. We generated particle number distribution 

functions for fugacities in the range z/zq G [0,1] by 
reweighting. 

Bimodal particle number distribution functions show 
a low-density peak that is given by the binomial distri- 
bution of non-interacting particles 

^-(-^= (l + .2U)^ (n)^^^^"'"^" ^''^ 
Thus, we identify the low-density peak as 

H^{n) ^ P^{n)^^ (59) 

where the factor P(0)/Pl(0) guarantees that the distri- 
bution function P{n) and the low-density peak H\,{n) co- 
incide for the empty tube, i.e., n = O. The second peak is 
a high-density, peak defined by Hii{n) = P{n) — Hi^{n). 

Next, we define a measure M that quantifies the ex- 
tent of bistability. It should be large if the areas be- 
low the two peaks, given by A^h = Sn-^H("-) and 
N-L — J2n H-L{n) = 1 — A^H, are of comparable size. Thus, 
we form the product of the two areas 

A = 4A^hAl (60) 

The factor 4 ensures that ^ = 1 if the two peaks have 
equal weight. If one peak dominates the particle number 
distribution function then A ^ 0. 

If the two peaks were due to coexistence at a first or- 
der phase transition, changing the fugacity would only 
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FIG. 14: Particle distribution functions for = 500 and 
zjzQ = 0.085 for the full Hamiltonian (solid, black) and in 
the SCA (solid, red). The dashed lines show fits of the low 
density peaks with binomial distributions of non-interacting 
single dipoles. The dashed-dotted line is a fit of a binomial 
distribution of non-interacting chains (with charges at their 
ends) up to a length of ten particles. The inset shows the 
distribution functions for the full range of the particle number. 



change the relative weight of the two peaks. Here, the 
location of the peaks changes with the fugacity. Thus, 
our measure should also include the distance between 
the two peaks. The positions of the peaks are given by 
the mean values of the low-density and the high-density 
peak and their distance R can be written as 



R 



N ^ 



(61) 



Thus, we define our measure for bistability AI as the 
product of A and the distance R 



M = RA = ARNhNi^ 



(62) 



Figure [12] shows a density plot of this measure as a 
function of the inverse system size and the chemical po- 
tential. For small systems, the bimodal structure of the 
particle number distribution function can be seen for am- 
bient conditions corresponding to n = jiq, as was also 
observed in computer simulations [14| . For larger system 
sizes the range of the chemical potential where bimodal- 
ity is observed becomes narrower and the bimodality it- 
self weaker. For N ^ 1000 bimodality vanishes com- 
pletely. Also shown are the lines where the system is 
empty (average density n/N = p = 0.01), full [p = 0.99), 
and half filled [p = 0.5). The maximum of the bistabil- 
ity measure is always to the left, i.e., at lower chemical 
potential, of the line of half filling. The bistability map 
for NCA (not shown) is nearly identical to the map for 
SCA. 

The situation for the full Hamiltonian is slightly dif- 
ferent to that for the SCA as exemplified in Fig.[T4l The 
low-density peak is not perfectly reproduced by a bino- 
mial distribution of non-interacting, single dipoles. In- 
stead, we have to take into account the contributions of 



longer chains, that do not interact with each other, but 
whose energy is lowered (compared to the SCA) due to 
the effective charges at their ends. The peak is well re- 
produced if we include chains with lengths up to ten sites. 
Thus, the form of the low density peak can be explained 
by the lower energies of chains in the full Hamiltonian 
compared to the SCA. 



VI. SUMMARY AND CONCLUSIONS 

We have studied a onc-dimensional dipole model in 
various equivalent representations, and also in approx- 
imate formulations. The dipole model is a simple and 
powerful description of water in nanopores. It allows us 
to characterize the phase-like properties of confined wa- 
ter up to macroscopic dimensions [l7j . Here we have 
explored the filling thermodynamics by calculating ad- 
sorption isotherms and particle number fluctuations as a 
function of relative humidity. We have also characterized 
the regime of bistability, with distinct gas-likc and liquid- 
like states, as a function of tube length and fugacity. 

The model can be represented in terms of individ- 
ual dipoles, ordered segments, or effective charges with 
Coulomb-like interactions. The segment picture is not 
only an essential step in the derivation of the charge rep- 
resentation from the dipole picture; it also results in a 
simple characterization of the configuration space and 
thus an efficient formulation of non-local trial moves for 
Monte Carlo simulations. The charge picture is the phys- 
ically most appealing representation as the long-range 
interactions are due to charges at the ends of chain seg- 
ments. It also reduces the computational cost of calcu- 
lating the Hamiltonian. 

In the charge picture, the Coulomb-like effective in- 
teractions of defects and chain ends result from effective 
charges, whose magnitude depend on the dipole moment 
of water molecules and their distance in the water wire. 
L- and D- defects carry effective charges of opposite sign 
with a magnitude that is twice that of the charges at the 
chain cndpoints. An excess proton in the chain, corre- 
sponding to an L-dcfcct with an extra proton, carries an 
effective charge that is reduced considerably with respect 
to the charge of the bare proton. This polarization ef- 
fect is important for the free energetics of proton transfer 
through water-filled narrow pores Q . 

We introduced two approximate representations of the 
dipole model. By neglecting the long-range interaction in 
the charge representation, one arrives at a Hamiltonian 
that only depends on the particle number, the chain num- 
ber, and the defect number. This approximation can be 
improved by treating charge interactions within segments 
of length one explicitly. The resulting approximation pro- 
duces accurate descriptions for observables depending on 
the numbers of particles and fragments. For this so-called 
singlet-charge-approximation we can count the number of 
states as a function of the independent variables. This 
model can thus be formulated in terms of an effective 
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Hamiltonian with a phase-space of only four dimensions. 

We used these approximations to study the bistabihty 
of the particle number distribution function, which gets 
weaker with increasing system size. The biggest effect of 
the long-range Coulomb-like interactions, as compared to 
the short-range interactions in an Ising-like lattice gas, is 
to lower the energies of short chains. This energy low- 
ering accounts for the low-density peaks in the particle- 
number distributions, with small but non-zero particle 
numbers. Nevertheless, for system sizes that are not too 
small, the bistable behavior of water in nanopores is cap- 
tured by the SCA model without long-range interactions. 

The dipole model introduced here is general and should 
be applicable to other quasi-one dimensional systems 
with dipolar interactions, including polar fluids other 
than water as well as magnetic nanoparticlcs and colloids 
[H, [13 ■ It should also prove useful in studies of threc- 
dimensionally packed arrays of one-dimensional chains, 
such as those formed in membranes of parallel nanochan- 
nels 
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APPENDIX A: NUMBER OF STATES 

For the NCA and the SCA we can count the number 
of states as a function of the independent variables of 
the respective Hamiltonians. Here, we show the deriva- 
tion for the SCA, where the number of states depends on 
the particle number n, chain number ric, defect number 
rid, and number ni of single dipoles corresponding to seg- 
ments of length one. For the simpler NCA, the derivation 
of the number of states is similar. 

To model the free boundary conditions we add to our 
fully occupied system consisting of N sites an empty sec- 
tion of length one at each end. Thus, the fully occupied 
system is given by an empty section of length one, an 
occupied section of length iV, and again an empty sec- 
tion of length one. These end sections are useful for the 
formulation of trial moves and for the derivation of the 
number of states, as we will see in the following. 

To calculate the number of states we first take a section 
of length n — Hd — ni and split it into rig — ni parts, each 



of which is at least of length two, i.e., Zmin = 2. To split a 
section of M sites into m parts we have to choose m — 1 of 
the M — 1 points between sites where the section can be 
split. The number of possibilities to do so is given by the 
binomial coefficient (^^ l) . If we demand that each new 

\m—l J 

section consists of at least l^un sites then the number of 
points where we can split the section into parts is reduced 
by (Zmin — l)m. Thus, the number of possibilities, 7, of 
splitting a section with M sites in m sub-sections, where 
each section has a length I > Zmin is given by 



7 



M 



{Imin - i)m 
TO — 1 



1 



(Al) 



Inserting M — n — — ni, m = Ug — ni, and ^min = 2 in 
the above equation we obtain the first binomial coefficient 
of Eq. (122]). 

Next we count in how many ways we can combine 
the Tig — ni segments with lengths larger than one and 
the ni segments of length one to a particular sequence 
of these Ug segments. We do so by choosing positions 
for the ni identical, single dipoles out of Ug possible po- 
sitions, which gives the second binomial coefficient of 
Eq. (j33p . Then, we put the remaining Us — ni segments of 
lengths larger than one on the remaining positions with- 
out changing their order, which is uniquely determined. 

Grouping these segments to chains corresponds to 
splitting this sequence of rig segments into Uc parts, each 
consisting of at least one segment. This gives the third 
binomial coefficient of Eq. ((33|) . Each of these ric chains 
has two possible orientations which gives the factor 2""=. 

Finally, N — n + 2 empty sites have to be partitioned 
in 71c + 1 sections that are at least of length one, i.e., 
^min = 1 in Eq. (|A1[) . which gives the last Binomial co- 
efficient of Eq. (|22|) . By alternating empty sections and 
chains a particular configuration is obtained, which is 
again uniquely determined. 

To do Monte Carlo simulation of the SCA effective 
Hamiltonian given by Eq. (j45p . we need not only the 
degeneracy, but we also have to know the limits of the 
volume spanned by the variables {n, rid, rii}, i.e., their 
minimum and maximum values. The minimum particle 
= and the maximum particle number 



IS n 



max 



(A^) = N. The minimum number of defects is 
n™'" — and the maximum number is given by rid ~ 
for n < 3 and otherwise 



1 for n even 
for n odd 



The minimum number of chains is 

''0 for n = 
1 for n > 



n-'"(n) 



(A2) 



(A3) 



and the maximum number is n'^^'^^^N, 0, n^) = for the 
empty system and otherwise 

n^'"'{N, n, rid) = min {A^ - n + 1, Ji - 2?id} . (A4) 
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The minimum number of chains of length one is given by 

miiv \ J for n - nd > 2ns 

rij (n, rid, nc) = < 

I 27is — n + for n — rid < 2ns • 

(A5) 

The maximum number of chains of length 1 is given by 

nr-(n,nd,n.) = | ^ n - n, ^ ^^^^ 

ns — 1 for < n . 



APPENDIX B: TRIAL MOVES 

In the following we present the trial moves for the 
Monte Carlo simulation of the dipole lattice model. 

In the Metropolis algorithm [2^ the transition prob- 
ability from an old state o to a new state n is given 
by the product of the generation probability of a move, 
Pgcn{o n) and the acceptance probability, Pacc(o — > 
n). Imposing detailed balance, the Metropolis acceptance 
probability in the canonical ensemble is given by 



Piicc{o —nn) = min <^ 1 



Pgc„(o ^ n)e-PE(°) , 



(Bl) 



and correspondingly for Pacc{n o). In simulations of 
the grand-canonical ensemble, the energies contain addi- 
tional terms — nlnz where n is the fluctuating particle 
number. Usually the generation probabilities of the for- 
ward and the backward move are chosen to be equal and 
they cancel each other in the above equation. 

In our simulation, a configuration is given in the seg- 
ment picture, i.e., by the lengths of all sections and their 
orientations. This is a simple way to include the configu- 
rational constraints mentioned above but has the disad- 
vantage that for some trial moves the generation prob- 
ability for the forward and the backward move is asym- 
metric. These asymmetric generation probabilities have 
to be explicitly included in Eq. (|B1[) . This is the case for 
defect generation and recombination, chain splitting and 
joining, and the insertion and removal of a single dipole, 
as is explained below. 

For simplicity, we use in the remaining part of this sec- 
tion the word "choose" when we mean "choose with equal 
probability", i.e., when we draw some quantity from a 
uniform distribution. 

For the displacement of a defect we choose a chain c £ 
{1, . . . ,nc} that consists of m segments from which we 
choose segment i G {l,...,rn — 1}. We change the length 
of segment i by dAl and the length of segment i + 1 by 
—dAl where we have chosen a length G {1, Amax} and 
a direction d = ±1. The generation probability is given 
by 



pdis _ 
gen 



1 



2nc(m - 1)A„ 



(B2) 



For the generation of a defect we choose a chain c G 
{1, . . . ,nc} that consists of m segments from which we 



choose segment i £ {1, . . . , m}. If the length of this seg- 
ment, /, is long enough to carry a defect, i.e., / > 3, then 
we choose a length I' G {1, . . . , ^ — 2} and a direction 
d = ±1. The segment i is split in two chains of length 
h = I' and I2 = 1 — 1' and all chains of fragment /, on the 
side given by d are reoriented. The generation probability 
is given by 



pgon _ 2 

2nc7n(; - 2) 



(B3) 



This move increases the number of defects nd by one. 

For defect recombination we choose a chain c G 
{1, . . . , nd that consists of m segments from which we 
choose segment i G {l,...,7n — 1}. Additionally we 
choose a direction d ~ ±1 and join segments i and i-f 1 to 
a new segment with length I = li + /i+i -I- 1. All segments 
on the side d are reoriented. The generation probability 
is given by 



prcc 



1 



2nc{m — 1) 



(B4) 



This moves decreases the number of defects nd by one. 

For the displacement of a fragment we choose a chain 
c G {l,...,nc}, a direction = ±1, and a displacement 
A/ G {1, . . . , Amax}- The empty section on the side d of 
the fragment is lengthened by A^ and the empty section 
on the opposite side is shortened by AL The generation 
probability is given by 



pira. _ 



1 



^"^^ 2n,A, 



(B5) 



The generation probability for the reorientation of a 
chain, i.e., the reorientation of all segments of chain c G 
{1, . . . ,nc}, is given by 



prco 
gen 



(B6) 



The exchange move shortens a segment i G {1, . . . , mi} 
of fragment ci consisting of mi segments, and lengthens a 
segment j G {1, . . . , 7712} of C2 consisting of m2 segments a 
length A^ G {1, Amax}, therefore conserving the number 
of occupied sites. The generation probability is given by 



pcxc 

gon — 



1 



(B7) 



To split a chain in two chains we choose a chain ci G 
{1, . . . , nc — 1} and a segment j G {1, . . . , m} of length /. 
Next we have to choose where in the segment i a bond is 
broken by choosing a length V G {1, Z— 1}. One of the new 
fragments is displaced (d = ±1) a length AZ G {1, A,„ax}- 
The generation probability is given by 



pspi _ \, 

2n^m{l - 1)A„ 



(B8) 



This move increases the number of chains by one and the 
number of segments by one. 
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The inverse move is the joining of two chains to a single 
chain. We choose a chain c £ {1, . . . , 7t,c — 1}. If the last 
ordered segment of chain c and the first ordered segment 
of chain c + 1 have the same direction and if they are not 
further apart than A^ax (i-e., the length of the empty 
section between them is I < Amax) then we try to join 
them. We choose a direction c? = ±1 that decides if the 
left chain is moved towards the right or the right towards 
the left. Wc get for the generation probability 



pjoi _ 
gen 



1 



2(nc - 1) 



(B9) 



This move decreases the number of chains by one and the 
number of ordered segments by one. 

Next we present moves that change the occupation 
number. The transfer move adds or removes dipoles at 
the end of chains. First we choose a chain cG {l,...,nc}, 
at which end particles are transfered {di = ±1). and how 
many particles {Al G {1, Amax}) are either added or re- 
moved by lengthening or shortening of the chosen end 
segment {d2 = ±1)- Lengthening of the end segment 
is only possible if the empty section next to it is longer 
than the number of added particles. This also guarantees 
that the number of chains is not changed. This gives a 
generation probability of 



ptra |_ 



(BIO) 



The above move is only applicable if there are already 
occupied sites. Therefore, we also insert single dipoles 
in empty sections. To do so we choose an empty section 
i G {l,...,nc-l-l} and a site by choosing a length I' G 
{1,^ — 2} for the empty section on the left of the inserted 
dipole, which we assign an orientation d = ±1. This 
results in a generation probability 



pins 
gen 



1 



2(nc + !)(;- 2) 



(Bll) 



This move increases the occupation number, the number 
of chains, and the number of ordered segments by one. 

The inverse move removes a single dipole by choosing 
a chain c G {1, . . . ,nc} and checking if it is of length 
one. If so, we remove the single dipole by eliminating 
this chain consisting of segment j and the empty section 
to its right with index j + 1 . The empty section to its left 
gets the new length lj_i = lj^i+lj+i + 1. The generation 
probability is 



tribution function Pw{n), that stems from a biased sam- 
pling simulation at the fugacity z, by unfolding of the 
weight function w(n) and reweighting to the new fugac- 
ity z', 



P(n,/)==^P4n).--^"^ 
with the normalization constant N{z') given by 



N{z') 



N 

E 

n=0 



P„(n)z-"'(") 



(CI) 



(C2) 



If we want to calculate the average value of the observ- 
able O (which can be any element of the sampled list or 
a function of these elements) from a biased simulation, 
we need the joint distribution function of the order pa- 
rameter and the observable in the biased ensemble given 
by 

Pw {n, 0) = ^J2 '^("^'■' - n)S{0'-'^ - O) (C3) 

i 

for a discrete observable O, where Al is the number of 
samples and S{x) is Dirac's delta function. Wc obtain the 
average of the observable at a fugacity z' by evaluating 

, (C4) 



<'5>-]4i:i:^-(">'5)'5^-'"'(f)' 



Instead of calculating the two-dimensional histogram 
and performing the above average, we calculate the fol- 
lowing average of the observable O for each value of the 
order parameter in the biased ensemble, i.e.. 



(0)u:(n) 



o 



(C5) 

We then do the unfolding of the weight function and 
the reweighting to the new fugacity z' in a single step and 
obtain for the average of the observable O as a function 
of the order parameter 

The average value of the observable O at the new fugacity 
z' is then obtained as {O) = X]n(^)('^)- 



APPENDIX D: LATTICE GAS 



prcm 



(B12) 



This move decreases the occupation number, the number 
of chains, and the number of ordered segments by one. 



APPENDIX C: BIASED SAMPLING 

The particle number distribution function P(n, z') at 
the fugacity z' is obtained from the particle number dis- 



We determine the chemical potential, /Ui/2, where the 
system is half filled by exploiting the isomorphism be- 
tween the Ising model in an external field in the canon- 
ical ensemble and a lattice gas in the grand canonical 
ensemble. This allows us to determine ^1/2 from symme- 
try considerations. The Hamiltonian of the lattice gas is 
given by 



H = ^^ninjJ(|j -i|) 



(DI) 
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where = 0, 1 is the occupation number of site i and 
the interaction potential given by 



for i = j 

for \3-i\^l (D2) 
else . 



The grand canonical partition function of this lattice gas 
is given by 



S= ^ exp 

)ifc=0,l 



-/3 (^H-^iY.^}j 



(D3) 



which is isomorphic to the canonical partition function 
of the Ising model with dipole-dipole interactions 



Q 



with Si = 1 corresponding to = 1, = — 1 correspond- 
ing to ?ij = 0, J{\j - i\) = J{\j - z|)/4, and ^ = 4 J + 2ft. 

with J = X^fc^i '^('i^)- For vanishing external field h the 
magnetization of the Ising model vanishes which corre- 
sponds to a half filled state for the lattice gas. Thus, 
Ml/2 = 4 J which gives 



Ml/2 = E, + [1- C(3)] 



(D5) 



(D4) in agreement with Eq. ([SO]) . 
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